## Agroforestry Paper Analysis import
## tree dataset
library(readxl)
nyando1 <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "set")
View(nyando1)
# Create for Total Lower Project
nyandolp <- subset(nyando1, GroupTypeInRegion == 
    "Lower Project")
View(nyandolp)

# Create for Total Lower Control
nyandolc <- subset(nyando1, GroupTypeInRegion == 
    "Lower Control")
View(nyandolc)

# Create for Total Middle Project
nyandomp <- subset(nyando1, GroupTypeInRegion == 
    "Middle Project")
View(nyandomp)

# Create for Total Middle Control
nyandomc <- subset(nyando1, GroupTypeInRegion == 
    "Middle Control")
View(nyandomc)

## Relate with Excel File..Count
## number of species Lower Project
tapply(nyandolp$Species, nyandolp$`Agroforestry Practice`, 
    summary)

cit1 <- nyandolp %>% group_by(`Agroforestry Practice`, 
    Species) %>% summarize(n())

## Lower Control
cit2 <- nyandolc %>% group_by(`Agroforestry Practice`, 
    Species) %>% summarize(n())
write.csv(cit2, "cit2.csv")

## Middle Project
cit3 <- nyandomp %>% group_by(`Agroforestry Practice`, 
    Species) %>% summarize(n())
write.csv(cit3, "cit3.csv")

## Middle Control
cit4 <- nyandomc %>% group_by(`Agroforestry Practice`, 
    Species) %>% summarize(n())
write.csv(cit4, "cit4.csv")

## Import data for Shannon index
## calculations Boundary Planting
library(readxl)
boundary <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "boundary")
View(boundary)

# Hedge Rows
hedgerow <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "hedgerow")
View(hedgerow)

# Riparian Buffers
riparian <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "riparian")
View(riparian)

# Woodlots
woodlot <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "woodlot")
View(woodlot)

## MPTs
mpt <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "mpt")
View(mpt)

################################################################# 

### Calculate Shannon Index
library(vegan)

# Boundary Planting ##1 Lower
# Project
lpbou <- boundary[c(1:25), c(1:2)]
lpbou$h <- lpbou$n1/sum(lpbou$n1)
lpbou$hp <- lpbou$h * log(lpbou$h)
-sum(lpbou$hp)
lpbou = t(lpbou)
lpbou = as.data.frame(lpbou)
colnames(lpbou) <- as.character(unlist(lpbou[1, 
    ]))
lpbou = lpbou[-1, ]
# Remove row names
rownames(lpbou) <- c()

lpbou = as.numeric(lpbou)

diversity(lpbou, index = "shannon")

## Lower Control
lcbou <- boundary[c(1:33), c(3:4)]
lcbou$h <- lcbou$n2/sum(lcbou$n2)
lcbou$hp <- lcbou$h * log(lcbou$h)
-sum(lcbou$hp)
lcbou = t(lcbou)
lcbou = as.data.frame(lcbou)
colnames(lcbou) <- as.character(unlist(lcbou[1, 
    ]))
lcbou = lcbou[-1, ]
# Remove row names
rownames(lcbou) <- c()

lcbou = as.numeric(lcbou)

diversity(lcbou, index = "shannon")

## Middle Project
mpbou <- boundary[c(1:51), c(5:6)]
mpbou$h <- mpbou$n3/sum(mpbou$n3)
mpbou$hp <- mpbou$h * log(mpbou$h)
-sum(mpbou$hp)
mpbou = t(mpbou)
mpbou = as.data.frame(mpbou)
colnames(mpbou) <- as.character(unlist(mpbou[1, 
    ]))
mpbou = mpbou[-1, ]
# Remove row names
rownames(mpbou) <- c()

mpbou = as.numeric(mpbou)

diversity(mpbou, index = "shannon")

## Middle Control
mcbou <- boundary[c(1:39), c(7:8)]
mcbou$h <- mcbou$n4/sum(mcbou$n4)
mcbou$hp <- mcbou$h * log(mcbou$h)
-sum(mcbou$hp)
mcbou = t(mcbou)
mcbou = as.data.frame(mcbou)
colnames(mcbou) <- as.character(unlist(mcbou[1, 
    ]))
mcbou = mcbou[-1, ]
# Remove row names
rownames(mcbou) <- c()

mcbou = as.numeric(mcbou)

diversity(mcbou, index = "shannon")



# Hedge Rows ####2 Lower Project
lphed <- hedgerow[c(1:16), c(1:2)]
lphed$h <- lphed$n1/sum(lphed$n1)
lphed$hp <- lphed$h * log(lphed$h)
-sum(lphed$hp)
lphed = t(lphed)
lphed = as.data.frame(lphed)
colnames(lphed) <- as.character(unlist(lphed[1, 
    ]))
lphed = lphed[-1, ]
# Remove row names
rownames(lphed) <- c()

lphed = as.numeric(lphed)

diversity(lphed, index = "shannon")

## Lower Control is zero

## Middle Project
mphed <- hedgerow[c(1:35), c(5:6)]
mphed$h <- mphed$n3/sum(mphed$n3)
mphed$hp <- mphed$h * log(mphed$h)
-sum(mphed$hp)
mphed = t(mphed)
mphed = as.data.frame(mphed)
colnames(mphed) <- as.character(unlist(mphed[1, 
    ]))
mphed = mphed[-1, ]
# Remove row names
rownames(mphed) <- c()

mphed = as.numeric(mphed)

diversity(mphed, index = "shannon")

## Middle Control is zero

## Woodlots #####3 Lower Project
lpwod <- woodlot[c(1:24), c(1:2)]
lpwod$h <- lpwod$n1/sum(lpwod$n1)
lpwod$hp <- lpwod$h * log(lpwod$h)
-sum(lpwod$hp)
lpwod = t(lpwod)
lpwod = as.data.frame(lpwod)
colnames(lpwod) <- as.character(unlist(lpwod[1, 
    ]))
lpwod = lpwod[-1, ]
# Remove row names
rownames(lpwod) <- c()

lpwod = as.numeric(lpwod)

diversity(lpwod, index = "shannon")

## Lower Control

lcwod <- woodlot[c(1:2), c(3:4)]
lcwod$h <- lcwod$n2/sum(lcwod$n2)
lcwod$hp <- lcwod$h * log(lcwod$h)
-sum(lcwod$hp)
lcwod = t(lcwod)
lcwod = as.data.frame(lcwod)
colnames(lcwod) <- as.character(unlist(lcwod[1, 
    ]))
lcwod = lcwod[-1, ]
# Remove row names
rownames(lcwod) <- c()

lcwod = as.numeric(lcwod)

diversity(lcwod, index = "shannon")

## Middle Project
mpwod <- woodlot[c(1:36), c(5:6)]
mpwod$h <- mpwod$n3/sum(mpwod$n3)
mpwod$hp <- mpwod$h * log(mpwod$h)
-sum(mpwod$hp)
mpwod = t(mpwod)
mpwod = as.data.frame(mpwod)
colnames(mpwod) <- as.character(unlist(mpwod[1, 
    ]))
mpwod = mpwod[-1, ]
# Remove row names
rownames(mpwod) <- c()

mpwod = as.numeric(mpwod)

diversity(mpwod, index = "shannon")

## Middle Control is zero

## MPTs (Multi Purpose) ###4 Lower
## Project
lpmpt <- mpt[c(1:30), c(1:2)]
lpmpt$h <- lpmpt$n1/sum(lpmpt$n1)
lpmpt$hp <- lpmpt$h * log(lpmpt$h)
-sum(lpmpt$hp)
lpmpt = t(lpmpt)
lpmpt = as.data.frame(lpmpt)
colnames(lpmpt) <- as.character(unlist(lpmpt[1, 
    ]))
lpmpt = lpmpt[-1, ]
# Remove row names
rownames(lpmpt) <- c()

lpmpt = as.numeric(lpmpt)

diversity(lpmpt, index = "shannon")

## Lower Control
lcmpt <- mpt[c(1:31), c(3:4)]
lcmpt$h <- lcmpt$n2/sum(lcmpt$n2)
lcmpt$hp <- lcmpt$h * log(lcmpt$h)
-sum(lcmpt$hp)
lcmpt = t(lcmpt)
lcmpt = as.data.frame(lcmpt)
colnames(lcmpt) <- as.character(unlist(lcmpt[1, 
    ]))
lcmpt = lcmpt[-1, ]
# Remove row names
rownames(lcmpt) <- c()

lcmpt = as.numeric(lcmpt)

diversity(lcmpt, index = "shannon")

## Middle Project
mpmpt <- mpt[c(1:67), c(5:6)]
mpmpt$h <- mpmpt$n3/sum(mpmpt$n3)
mpmpt$hp <- mpmpt$h * log(mpmpt$h)
-sum(mpmpt$hp)
mpmpt = t(mpmpt)
mpmpt = as.data.frame(mpmpt)
colnames(mpmpt) <- as.character(unlist(mpmpt[1, 
    ]))
mpmpt = mpmpt[-1, ]
# Remove row names
rownames(mpmpt) <- c()

mpmpt = as.numeric(mpmpt)

diversity(mpmpt, index = "shannon")

## Middle Control
mcmpt <- mpt[c(1:58), c(7:8)]
mcmpt$h <- mcmpt$n4/sum(mcmpt$n4)
mcmpt$hp <- mcmpt$h * log(mcmpt$h)
-sum(mcmpt$hp)
mcmpt = t(mcmpt)
mcmpt = as.data.frame(mcmpt)
colnames(mcmpt) <- as.character(unlist(mcmpt[1, 
    ]))
mcmpt = mcmpt[-1, ]
# Remove row names
rownames(mcmpt) <- c()

mcmpt = as.numeric(mcmpt)

diversity(mcmpt, index = "shannon")

# Riparian Buffers ####5 Lower
# Project
lprip <- riparian[c(1:7), c(1:2)]
lprip$h <- lprip$n1/sum(lprip$n1)
lprip$hp <- lprip$h * log(lprip$h)
-sum(lprip$hp)
lprip = t(lprip)
lprip = as.data.frame(lprip)
colnames(lprip) <- as.character(unlist(lprip[1, 
    ]))
lprip = lprip[-1, ]
# Remove row names
rownames(lprip) <- c()

lprip = as.numeric(lprip)

diversity(lprip, index = "shannon")

# Lower Control is zero Medium
# Project is zero

## Middle Control
mcrip <- riparian[c(1:2), c(7:8)]
mcrip$h <- mcrip$n4/sum(mcrip$n4)
mcrip$hp <- mcrip$h * log(mcrip$h)
-sum(mcrip$hp)
mcrip = t(mcrip)
mcrip = as.data.frame(mcrip)
colnames(mcrip) <- as.character(unlist(mcrip[1, 
    ]))
mcrip = mcrip[-1, ]
# Remove row names
rownames(mcrip) <- c()

mcrip = as.numeric(mcrip)

diversity(mcrip, index = "shannon")

################################################## 

## Calculating Shannon Indices For
## Sites and Whole region
library(readxl)
shannon <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "shannon")
View(shannon)

## Lower Site
lshan <- shannon[c(1:48), c(1:2)]
lshan$h <- lshan$n1/sum(lshan$n1)
lshan$hp <- lshan$h * log(lshan$h)
-sum(lshan$hp)
lshan = t(lshan)
lshan = as.data.frame(lshan)
colnames(lshan) <- as.character(unlist(lshan[1, 
    ]))
lshan = lshan[-1, ]
# Remove row names
rownames(lshan) <- c()

lshan = as.numeric(lshan)

diversity(lshan, index = "shannon")

## Middle Site
mshan <- shannon[c(1:95), c(3:4)]
mshan$h <- mshan$n2/sum(mshan$n2)
mshan$hp <- mshan$h * log(mshan$h)
-sum(mshan$hp)
mshan = t(mshan)
mshan = as.data.frame(mshan)
colnames(mshan) <- as.character(unlist(mshan[1, 
    ]))
mshan = mshan[-1, ]
# Remove row names
rownames(mshan) <- c()

mshan = as.numeric(mshan)

diversity(mshan, index = "shannon")

## Whole Region
rshan <- shannon[c(1:102), c(5:6)]
rshan$h <- rshan$n3/sum(rshan$n3)
rshan$hp <- rshan$h * log(rshan$h)
-sum(rshan$hp)
rshan = t(rshan)
rshan = as.data.frame(rshan)
colnames(rshan) <- as.character(unlist(rshan[1, 
    ]))
rshan = rshan[-1, ]
# Remove row names
rownames(rshan) <- c()

rshan = as.numeric(rshan)

diversity(rshan, index = "shannon")

## Testing Diverse Package
## Calculating Shannon Indices For
## Sites and Whole region
library(readxl)
diverse <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "diverse")
View(diverse)
### Convert from list
diversea <- do.call(rbind.data.frame, 
    diverse)
diversity(diversea)
#### Seta a Block-Practice-Value
diversea = diverse[, c(1, 3, 6)]
diversity(diversea)


## Calculating Shannon Indices For
## Site Types
library(readxl)
typeshan <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "typeshan")
View(typeshan)
## Lower Project
lptshan <- typeshan[c(1:41), c(1:2)]
lptshan$h <- lptshan$n1/sum(lptshan$n1)
lptshan$hp <- lptshan$h * log(lptshan$h)
-sum(lptshan$hp)

## Lower Control
lctshan <- typeshan[c(1:38), c(3:4)]
lctshan$h <- lctshan$n2/sum(lctshan$n2)
lctshan$hp <- lctshan$h * log(lctshan$h)
-sum(lctshan$hp)

## Middle Project
mptshan <- typeshan[c(1:84), c(5:6)]
mptshan$h <- mptshan$n3/sum(mptshan$n3)
mptshan$hp <- mptshan$h * log(mptshan$h)
-sum(mptshan$hp)

## Middle Control
mctshan <- typeshan[c(1:67), c(7:8)]
mctshan$h <- mctshan$n4/sum(mctshan$n4)
mctshan$hp <- mctshan$h * log(mctshan$h)
-sum(mctshan$hp)

## Calculating Community Dataset
library(readxl)
anovan <- read_excel("C:/Users/LORERO/OneDrive - CGIAR/Desktop/2019 Nyando Data.xlsx", 
    sheet = "anova")
View(anovan)
# Extract dataset
anovana = anovan[c(1:4), c(18:119)]
View(anovana)
rownames(anovana) = c("lp", "lc", "mp", 
    "mc")
View(anovana)

diversityresult(anovana, y = NULL, index = c("Shannon"), 
    method = "each site", sortit = TRUE, 
    digits = 5)
diversityresult(anovana, y = NULL, index = c("richness"), 
    method = c("each site", "pooled"), 
    sortit = TRUE, digits = 5)
diversityresult(anovana, y = NULL, index = c("abundance"), 
    method = "each site", sortit = TRUE, 
    digits = 5)
diversityresult(anovana, y = NULL, index = c("Simpson"), 
    method = "each site", sortit = TRUE, 
    digits = 5)
diversityresult(anovana, y = NULL, index = c("Eevenness"), 
    method = "each site", sortit = TRUE, 
    digits = 5)



nyando_agb <- nyando1
nyando_agb$number <- rep(1:length(nyando1$IntervieweeName))#add count
nyando_agb_main <- nyando_agb[, c("Block", 
                                  "type_r", "Agroforestry_Practice", 
                                  "AGB (Mg)", "number")]
#aggregate dataset
aggdata_agb <- aggregate(. ~ number + Block + 
                           type_r + Agroforestry_Practice, 
                         nyando_agb_main, sum)

aggdata_agb$land <- nyando.socio$land_owned_ha[
  match(aggdata_agb$IntervieweeName, nyando.socio$name)]

aggdata_agb$agb = aggdata_agb$`AGB (Mg)`
aggdata_agb$agfp = aggdata_agb$`Agroforestry Practice`
aggdata_agb$Block = as.factor(aggdata_agb$Block)
aggdata_agb$type_r = as.factor(aggdata_agb$type_r)

## Lower vs Middle ANOVA for AGB
oneway.test(agb ~ Block, data = aggdata_agb)  #ANOVA Unequal var

## Lower Project AGB
aggdata_agb1 <- aggdata_agb %>% filter(type_r == "Lower Project")

# Total AGB under Boundary Planting
sum(filter(aggdata_agb1, agfp == "Boundary planting")$agb)


# Total AGB under Riparian buffers
sum(filter(aggdata_agb1, agfp == "Riparian buffers")$agb)
## [1] 0.09491625

## Lower Control AGB
aggdata_gb2 <- filter(aggdata_gb, type_r == 
                        "Lower Control")

## Middle Project AGB
aggdata_gb3 <- filter(aggdata_gb, type_r == 
                        "Middle Project")

# Middle Control AGB
aggdata_gb4 <- filter(aggdata_gb, type_r == 
                        "Middle Control")
